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ABSTRACT 

We present the SuperNova Explosion Code (SNEC), an open-source Lagrangian code for the hydrodynamics 
and equilibrium-diffusion radiation transport in the expanding envelopes of supernovae. Given a model of a 
progenitor star, an explosion energy, and an amount and distribution of radioactive nickel, SNEC generates 
the bolometric light curve, as well as the light curves in different broad bands assuming black body emission. 

As a first application of SNEC, we consider the explosions of a grid of 15 Mq (at zero-age main sequence) 
stars whose hydrogen envelopes are stripped to different extents and at different points in their evolution. The 
resulting light curves exhibit plateaus with durations of ~20 — 100 days if > 1.5 — 2Mq of hydrogen-rich 
material is left and no plateau if less hydrogen-rich material is left. If these shorter plateau lengths are not 
seen for Type IIP supernovae in nature, it suggests that, at least for zero-age main sequence masses < 20 Mq, 
hydrogen mass loss occurs as an all or nothing process. This perhaps points to the important role binary 
interactions play in generating the observed mass-stripped supernovae (i.e., Type Ib/c events). These light 
curves are also unlike what is typically seen for Type IIL supemovae, arguing that simply varying the amount 
of mass loss cannot explain these events. The most stripped models begin to show double-peaked light curves 
similar to what is often seen for Type lib supernovae, confirming previous work that these supernovae can 
come from progenitors that have a small amount of hydrogen and a radius of ~ 500 Rq. 

Subject headings: hydrodynamics — radiative transfer — supernovae: general 


1. INTRODUCTION 

During the last half century the number of observed super¬ 
novae (SNe) has increased exponentially (Minkowski 1964; 
Cappellaro 2014). Much of this progress has been fueled 
by recent surveys, such as the Lick Observatory Supernova 
Search (LOSS, Leaman et al. 2011; Smith et al. 2011), the 
Sloan Digital Sky Survey (SDSS, Frieman et al. 2008), and 
the Palomar Transient Factory (PTE, Rau et al. 2009). In addi¬ 
tion to providing more complete and detailed samples of well- 
known classes of SNe (Type la, Ib/c, II), these surveys have 
found a wide range of previously unknown explosive events, 
from superluminous SNe (Quimby et al. 2011) to rapid SN- 
like transients (Perets et al. 2010; Kasliwal et al. 2010, 2012; 
Foley et al. 2013; Inserra et al. 2015). This has opened our 
eyes to the broader range of astrophysical explosions that can 
exist in nature. 

Progress in explosive, transient observations has been 
closely followed by progress in analytic and numerical 
light curve modeling. For example, for SNe IIP, this has 
ranged from analytic scalings (Arnett 1980; Chugai 1991; 
Popov 1993) to detailed numerical works (e.g., Litvinova & 
Nadezhin 1983; Chieffi et al. 2003; Young 2004; Kasen & 
Woosley 2009; Bersten et al. 2011; Dessart et al. 2013). These 
investigations focused on understanding the general imprint 
of progenitor characteristics (mass, radius, abundance and 
mixing of ^®Ni, etc.) on the shape and luminosity of SN light 
curves. In other cases, detailed comparisons have been made 
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between specific SNe II and numerical models (Arnett 1988; 
Shigeyama et al. 1988; Woosley 1988; Utrobin 1993 for SN 
1987A, Nomoto et al. 1993; Bartunov et al. 1994; Shigeyama 
et al. 1994; Young et al. 1995; Blinnikov et al. 1998 for SN 
1993J, Baklanov et al. 2005; Utrobin 2007 for SN 1999em). 

The combination of growing samples of SNe and other pre¬ 
viously unknown transients has motivated us to develop a 
new code for numerical studies of the light curves of SN and 
SN-like explosions. Called the SuperNova Explosion Code 
(SNEC), this general purpose code will allow the user to take 
a stellar model (or other ad-hoc density profile with other 
thermodynamic and compositional information), input energy 
to generate an explosion, follow the hydrodynamic response, 
and produce light curves. The current iteration of SNEC is 
spherically-symmetric (ID), and uses Lagrangian hydrody¬ 
namics and equilibrium-diffusion (one-temperature) radiation 
transport. It also follows other basic physics needed for light 
curves such as ionization and heating from ®®Ni. 

In terms of complexity and amount of physics included, 
SNEC is at a somewhat intermediate position compared with 
existing SN light curve codes. The current state-of-the-art are 
multi-group radiation-hydrodynamic codes (as in Blinnikov 
& Bartunov 1993; Moriya et al. 2011) and line-transfer ra¬ 
diative transfer codes that assume homologous expansion and 
either make the local thermodynamic equlibrium approxima¬ 
tion (LTE; e.g., Kasen & Woosley 2009) or are fully non-LTE 
(e.g., Hillier & Dessart 2012). SNEC bridges the gap between 
these codes and analytical investigations, e.g., those of Arnett 
(1980), Chugai (1991), and Popov (1993), and the more re¬ 
cent ones of Nakar & Sari (2010), Goldfriend et al. (2014), 
and Nakar et al. (2015). Our work is very much in the same 
spirit as the works of Bersten et al. (2011) and Ergon et al. 
(2015). A crucial aspect of SNEC is that unlike these other 


2accepted to ApJ, 10/15/2015 


Morozova et al. 


codes, it is open source and publicly available"^. This will 
make light curve modeling accessible and reproducible for 
the broader community. It can be used for a wide range of 
studies, from generating typical SN light curves as an edu¬ 
cational tool to making light curves for novel explosion sce¬ 
narios. Modeling explosions and light curves involves a wide 
range of physics and necessary approximations. Hence, mak¬ 
ing code available and results reproducible is crucial for the 
advancement of the field. 

An important strength of SNEC is the ability to systemati¬ 
cally and quickly explore changes in stellar properties to learn 
how they impact the resulting light curves. This is especially 
useful for investigating the underlying mechanisms behind 
the photometric diversity of SN II light curves, such as the 
SNe IIP (with nearly constant plateau luminosity for a period 
^ 100 days past maximum, and the most common type), SNe 
IIL (with linearly declining magnitude), and SNe Ilb (which 
show signatures of hydrogen present, but with a light curve 
generally more similar to SNe Ib)^. In particular, the con¬ 
nection between SNe IIP and SNe IIL has long been a point 
of contention in the SN community. Early on, it was sug¬ 
gested by Barbon et al. (1979) and corroborated by Blinnikov 
& Bartunov (1993) that the morphological differences might 
be explained by altering the envelope masses while keeping 
the explosion mechanism the same (however, see Swartz et al. 
1991 for an alternative picture). In this explanation, SNe IIL 
would simply have less hydrogen-rich envelopes than SNe IIP. 
Nevertheless, SNe IIL must still have appreciable hydrogen 
present, otherwise they would become SNe Ib/c instead. This 
suggests that there may exist a continuous range of hydro¬ 
gen mass stripping and thus a continuous range of events be¬ 
tween canonical SNe IIP and IIL. Furthermore, SNe lib have 
been inferred to have a small amount of hydrogen present 
(~ 0.01 — 0.1 Mq, Woosley et al. 1994; Bersten et al. 2012; 
Nakar & Piro 2014), and thus in principle with sufficient mass 
loss a transition should be seen all the way to SNe lib. The 
question is whether additional ingredients are needed beyond 
just increased mass loss to reproduce these features. 

Motivated by these questions, we investigate the mass-loss 
hypothesis for the origin of these SN II classes as a first ap¬ 
plication of SNEC. We use presupernova stellar models gen¬ 
erated with the MESA code (Paxton et al. 2011, 2013). Be¬ 
sides providing excellent control in generating models and 
investigating mass stripping, MESA has been used for other 
light curve studies (e.g., Dessart et al. 2013), which allows 
us to compare directly as a further check of SNEC. Although 
we find that varying levels of hydrogen mass stripping short¬ 
ens the plateau of the light curves, we conclude that simply 
varying the amount of mass loss alone cannot explain the full 
range of properties of SNe IIL. In the most mass stripped 
cases, we begin to see double-peaked light curves reminiscent 
of some SNe lib, suggesting that this transition occurs more 
naturally. Further work will be needed for a more complete 
investigation of SNe lib properties. 

In Section 2, we describe SNEC in detail. We follow, in 
Section 3, with our study of massive stars with varying levels 
of stripping. In Section 4, we present the resulting SN light 
curves. We conclude in Section 5 with a summary of our find¬ 
ings and a discussion of future work. In Appendices A and B, 
we compare SNEC with two other codes. 

^ http://stellaroollapse.org/SNEC 

^ We ignore SNe tin for our study because they show clear signs of inter¬ 
action, which is beyond the scope of this work. 


2. SNEC 

We describe SNEC with a focus on SN IIP light curve 
modeling. Although we anticipate that SNEC will continue 
to evolve and improve as it is utilized for new projects and 
more physics is added, the discussion below will provide 
some necessary background and summarize SNEC’s gen¬ 
eral features. A more detailed description that also includes 
the finite-difference form of the equations and implementa¬ 
tion details is available on the SNEC webpage, http:// 
stellarcollapse.org/SNEC. 

It is also important to compare SNEC with other SN light 
curve codes. Although below we focus on what is imple¬ 
mented in SNEC, in Appendices A and B we consider the 
work of Bersten et al. (2011), whose code has a similar level 
of complexity as SNEC, and Dessart et al. (2013), whose code 
performs a more detailed treatment of the radiative transfer. 
We find that both comparisons give satisfactory results with 
the main difference being the transition from the plateau to the 
®®Ni tail found by Dessart et al. (2013). This disagreement 
likely reflects an intrinsic difference between equilibrium- 
diffusion radiation-hydrodynamics codes, such as SNEC, and 
line-transfer radiative-transfer codes such as CMFGEN used in 
Dessart et al. (2013). 


2.1. ID Lagrangian LTE Radiation Hydrodynamics with 
Ionization 

Lagrangian hydrodynamics in spherical symmetry, supple¬ 
mented with a radiation diffusion term, written to the first or¬ 
der m v/c (see, e.g., Mihalas & Mihalas 1984; Mezzacappa 
& Bruenn 1993; Bersten 2010), results in a mass conserva¬ 
tion equation (continuity equation). 
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Here m = p{r')dr' is the mass coordinate, r is the 

radius, t is the time, p is the density, e is the specific internal 
energy, P is the pressure, v = dr/dt is the velocity of the 
matter, Q is the artificial viscosity, and G is the gravitational 
constant. In Equation (2), cni is the specific energy deposi¬ 
tion rate due to the radioactive decay of ®®Ni, which is equal 
to the time-dependent rate of energy release per gram of ra¬ 
dioactive nickel Crad, multiplied by the deposition function d, 
both defined below in Section 2.3 by Equations (14) and (15), 
respectively. The radiative luminosity L is 


L = 


Xac dT^ 
3k dm ’ 


(4) 


where a is the radiation constant, c is the speed of light, A 
is the flux-limiter and k is the Rosseland mean opacity. For 
capturing shocks, we use a simple von Neumann-Richtmyer 
artificial viscosity (Von Neumann & Richtmyer 1950) 


Q^fcpidv/di)^ 


if dv/dl < 0 , 
otherwise , 


(5) 
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where I is an integer grid coordinate, and C = 2. Following 
Bersten et al. (201 1), we use the flux limiter of Levermore & 
Pomraning (1981), 


where 


6 + 3i? 

6 + 3i? + > 


47rr^ 


dT^ 

dm 


(6) 

(7) 


We discretize in mass and time following the scheme 
of Mezzacappa & Bruenn (1993). The mass conservation 
and momentum conservation equations are updated time- 
explicitly. For Equation (2), we use a semi-implicit scheme 
with an adjustable parameter 6 in the discretization of the 
derivative dL/dm that can vary from fully explicit (0 = 0, 
only the luminosity from the previous time step is used in the 
scheme) to fully implicit {0 = 1, only the luminosity at the 
next time step is used in the scheme). The derivative dT^ jdm 
is linearized in 8T. We use 0 = 1/2 for all simulations pre¬ 
sented in this paper. Using an initial guess for the temperature 
at the next time step, we iteratively solve for ST, inverting 
a tridiagonal matrix each time, until the fractional change in 
temperature is less than a set tolerance (10“^ in the current 
version of the SNEC). We do not take the dependence of the 
opacity k on temperature into account in the implicit update 
and rather use the opacity from the previous time step when 
solving for the temperature at the next step. 

SNEC assumes local thermodynamical equilibrium (LTE), 
imposing the same temperature for radiation and matter. This 
assumption is not valid at shock breakout and during and af¬ 
ter the transition phase from optically thick to optically thin 
ejecta. In SNe IIP, it is reliable only during the plateau phase 
of the light curve (see the discussions in Blinnikov & Bar- 
tunov 1993; Bersten et al. 2011). However, the comparison 
performed by Bersten (2010) between her LTE code and a 
multi-group code suggests satisfactory agreement along the 
entire light curve. 

As a boundary condition for Equation (3), we adopt P = 0 
at the center of the boundary cell, so half a grid cell outside 
the star. Eor Equation (2), we assume that the luminosity at 
the surface is equal to the luminosity at the closest interior 
grid point, i.e., that the diffusive term, dL/dm, at the outer 
boundary is equal to zero. At the inner boundary, we take 
the velocity and the bolometric luminosity to be zero. In the 
modeling of core-collapse SN light curves, the inner boundary 
is typically not at m = 0 due to the presence of a neutron star 
(or a black hole), which is excluded from the grid. Setting the 
inner velocity to zero excludes any possibility for fallback of 
material onto the remnant in our models. 

To close the system of hydrodynamic equations, we em¬ 
ploy the analytic equation of state (EOS) given by Paczyhski 
(1983), hereafter the Paczyhski EOS. The Paczyhski EOS 
contains contributions from radiation, ions, and electrons, and 
takes into account electron degeneracy approximately. We re¬ 
peat some of our model calculations with the Helmholtz EOS 
(Timmes & Arnett 1999; Timmes & Swesty 2000), which in¬ 
cludes a (tabulated) complete electron EOS, to test the ap¬ 
proximations made in the Paczyhski EOS. We find that differ¬ 
ences between the Paczyhski and Helmholtz EOS have negli¬ 
gible influence on the resulting light curves. 

In order to account for recombination, we supplement the 
Paczyhski EOS with a routine that solves the Saha equations 
in the non-degenerate approximation as proposed in Zaghloul 


et al. (2000). The set of Saha equations, together with the con¬ 
dition of charge neutrality and number conservation of nuclei 
of a given chemical element (enumerated by index k), may be 
combined into a single transcendental equation for the aver¬ 
age charge of element k with atomic number as 
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i = 0,1,..., {Zk - 1) , 


where i is the number of the ionization state (0 corresponds 
to the neutral atom), Zk is the atomic number of element k, 
Uk is the number density of element k, gk^i is the statistical 
weight of the /-th ionization state of element k, Ik,i is the 
(positive) ionization energy for the ionization process i 
(i-fl), me is the electron rest mass, h is Planck’s constant and 
fce is Boltzmann’s constant. Equation (8) is solved iteratively 
at each call to the EOS, after which the ionization fractions 
ak,i are found as 



Although one can consider as many elements as necessary at 
the expense of computational time, for the present work we 
focus on the ionization of hydrogen and helium. The specific 
internal energy is calculated as 

^ — t^ion ^e\ Crad ^^ion j (H) 


where eion, CeU and Crad are the contributions from ions, elec¬ 
trons, and radiation, respectively, and 


Zfe / i-l 

^k ^ ^ I ^k,i ^ ^ 

i=l \ i=0 

is the ionization energy whose zero point for each element k 
is the neutral atom. 




2.2. Opacities 

The Rosseland mean opacity k is an essential input to 
our light curve models. In the high temperature regime 
(103.75 < T < 10® ^K), we use the OPAL Type II 

opacity tables (Iglesias & Rogers 1996) for solar metallic- 
ity (Zq = 0.02 here). These tables allow for an increase 
in the mass fractions of two chosen metals (in our case, car¬ 
bon and oxygen) by deducting an amount of helium to keep 
the sum of the mass fractions equal to unity. At low tem¬ 
peratures (lO^-^ K < T < 10^ ® K), we use the tables of 
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Ferguson et al. (2005). These tables are available for so¬ 
lar composition, but not for enhanced carbon and oxygen 
mass fractions compatible with the OPAL tables. In the over¬ 
lap region between the OPAL and the Ferguson et al. tables 
(IOS.t^k <T < lO^'^K), we give preference to the low- 
temperature opacity tables, because they take into account the 
contribution from molecular lines (this contribution is not in¬ 
cluded in the OPAL tables). For carbon and oxygen enhanced 
compositions, there are regions of low temperature and den¬ 
sity for which opacity values are not available. In these re¬ 
gions, the opacity is generally most sensitive to temperature, 
and thus we set the opacity to the nearest value that is avail¬ 
able at the same temperature. Most of the opacities set this 
way are below the opacity floor (see the discussion below), 
so that this deficiency in the tables has a small impact on the 
light curve evolution. At worst, it may affect the transition be¬ 
tween the plateau and the ^®Ni dominated part of SN IIP light 
curves when the photosphere first moves into carbon/oxygen- 
rich regions. 

The Rosseland mean opacity that we obtain from the OPAL 
and Ferguson et al. tables does not describe all possible 
sources of opacity needed for simulating SN light curves. 
As has been argued in previous works (see Karp et al. 1977; 
Young 2004; Blinnikov 1996; Bersten et al. 2011), the tabu¬ 
lated Rosseland mean opacity calculated for a static medium 
may underestimate the contribution of the line opacities in the 
rapidly expanding matter of the exploding star, plus it does 
not contain possible non-thermal ionization and excitation by 
gamma rays. Due to these missing effects, it is common prac¬ 
tice to use a so-called opacity floor, effectively imposing a 
minimum possible value for the opacity. Presently, there is no 
universally agreed-upon prescription for how to choose this 
opacity floor for a given composition and velocity. In previous 
work, different values of the opacity floor were chosen based 
either on simplified physical arguments (e.g., Shigeyama & 
Nomoto 1990) or based on comparisons with results obtained 
with multi-group or line-transfer codes (e.g., Bersten et al. 
2011). For SNe IIP, the values of the opacity floor for the 
hydrogen-rich envelope and the metal-rich core of the star, as 
well as the location and shape of the transition between the 
core/envelope opacities, can strongly influence the shape of 
the resulting light curve. Qualitatively (as shown in Bersten 
2010 and confirmed by our simulations), a lower value of the 
opacity floor in the envelope of the star increases the plateau 
luminosity and decreases the duration of the plateau, and vice 
versa. The luminosity of the plateau and its duration are im¬ 
portant observed photometric quantities that are used for sta¬ 
tistical studies of SNe IIP (as, e.g., in Anderson et al. 2014). 
Hence, it is important to keep the uncertainties in the opacity 
and their propagation into variations of the light curve in mind 
when comparing modeling results with observations. 

In the work of Bersten et al. (2011) on SN IIP light curves, 
the opacity floor was set to 0.01 cm^ for the “envelope” 
and 0.24 cm^ for the “core.” Since Bersten et al. (2011) 
are not specific in defining what constitutes the core and en¬ 
velope, and because this prescription introduces large opac¬ 
ity discontinuities, we take a different approach in SNEC. We 
choose the opacity floor to be linearly proportional to metal- 
licity Z at each grid point, and set it to 0.01 cm^ g“^ for solar 
composition (Z = 0.02) and to 0.24 cm^ g“^ for a pure metal 
composition (Z = 1)®. Note that we do not include the opac- 

* For the SNe IIP we study here, we do not find pure-metal regions in our 
models due to mixing that we impose during the explosion as described in 


ity floor in the calculation of the optical depth and position of 
the photosphere (as in Bersten et al. 2011). This is justified 
by the fact that the opacity floor is used to account for line ef¬ 
fects, which have minor influence on the shape of most of the 
continuum spectrum. However, we note that in the blue part 
of the spectrum (e.g., in the U and B bands), the continuum 
may be affected by the numerous lines of iron-group elements 
(see, e.g.. Figure 8 of Kasen & Woosley 2009). 

2.3. Radioactive ®®Ni and Bolometric Luminosity 

Radioactive ®®Ni in core-collapse SNe is synthesized by ex¬ 
plosive nuclear burning of intermediate-mass elements during 
the first seconds of the SN explosion in the inner regions of the 
star. It is mixed outward by hydrodynamic instabilities trig¬ 
gered by the shock’s propagation through the envelope (see, 
e.g., Kifonidis et al. 2003, 2006; Wongwathanarat et al. 2015). 
The gamma rays, emitted in the ®®Ni —>• ®®Co —)■ ®®Fe de¬ 
cay process, diffuse and thermalize, providing an additional 
source of energy eni in Equation (2). 

The present version of SNEC does not include a nuclear 
reaction network and the amount and distribution of synthe¬ 
sized ®®Ni is provided by the user. While this is a technical 
limitation that will be removed in future versions of SNEC, 
one should keep in mind that nucleosynthetic yields are sensi¬ 
tive to (1) how the explosion is launched, (2) where (in mass 
and spatial coordinate) it is launched, (3) to uncertainties in 
the structure and composition of the layers in which explo¬ 
sive burning occurs (e.g.. Young & Fryer 2007). Specifying 
the ®®Ni yield by hand removes these uncertainties from our 
models and has the added benefit of allowing the user com¬ 
plete control of radioactive heating, which can be useful for 
exploring how it impacts light curves. 

For the gamma rays released in ®®Ni and ®®Co decay, we 
follow the gray transfer approximation of Swartz et al. (1995), 
solving the transfer equation in the form 

dl' 

-r- = I' — , (13) 

dr 

where r is the optical depth along a given ray, Xni is the mass 
fraction of ®®Ni, /' = (Att and I is the energy- 
integrated intensity. The effective gamma-ray opacity is as¬ 
sumed to be purely absorptive and independent of energy, 
K.y = 0.06 Ye cm^ g“^, where Ye is the electron fraction. The 
time-dependent rate of energy release per gram of radioactive 
nickel, Crad, is equal to 

Crad = 3.9 X 10^° exp(—f/rNi) -1-6.78 x 10®[exp(—f/rco) 
- exp(-f/TNi)] ergg"^s"^ , (14) 

where tni and rco are the mean lifetimes of ®®Ni and ®®Co, 
equal to 8.8 and 113.6 days, respectively. The local heating 
rate in each grid point is equal to eni = where d is the 

deposition function 

d = ^ I'dixj, (15) 

where w is the solid angle. 

We do not take into account the energy from positrons re¬ 
leased in the radioactive decay of ®®Co (which occurs for 19% 
of the decays). The total kinetic energy of positrons per ®®Co 
decay is ~ 0.12 MeV versus ^ 3.61 MeV emitted via gamma 

Section 4.1. 
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rays (see Nadyozhin 1994). Therefore, neglecting this con¬ 
tribution constitutes an error of order 3 — 4% in the overall 
energetics of the ®®Ni decay (Swartz et al. 1995). 

Finally, we calculate the observed bolometric luminosity 
as suggested in Young (2004). It consists of two parts, the 
luminosity at the photosphere and the luminosity due to the 
absorption of gamma rays from ®®Ni/®®Co decay above the 
photosphere 

I'M 

^ohs — .^ptioto “f / Sfiep(^17l)dTTl . (16) 

•'mphoto 

Here TOphoto is the mass coordinate of the photosphere, M 
is the total mass of the star, S'dep is the energy per gram per 
second deposited by gamma rays. The location of the photo¬ 
sphere is defined by the optical depth r = 2/3, and Tphoto is 
found from Equation (4) at the photosphere location. 

3. PROGENITOR MODELS WITH VARYING 
HYDROGEN-RICH ENVELOPE MASSES 

3.1. Motivation and Overall Strategy 

As a first application of SNEC we investigate the light 
curves of SNe from massive stars that have lost varying 
amounts of their hydrogen-rich envelope during their evolu¬ 
tion. We use SNEC to explode presupernova models that we 
generate with the open-source stellar evolution code MESA 
(release version 6794; Paxton et al. 2011, 2013). The em¬ 
ployed MESA input files and final presupernova profiles are 
available at http : //stellarcollapse . org/SNEC. 

Massive stars may lose large fractions of their hydrogen- 
rich envelopes via steady line-driven winds, via stable or 
unstable binary mass transfer, or, possibly, through pulsa- 
tional instabilities and eruptions (see, e.g.. Smith 2014 for 
a recent review). The demographics of SN types combined 
with initial-mass-function considerations suggest that line- 
driven winds alone cannot account for the fraction of ob¬ 
served stripped-envelope SNe (Smith et al. 201 1). One of the 
other avenues of mass loss may be required to partially or 
completely remove hydrogen-rich envelope to account for the 
fraction of observed Type Ilb and Ib/c SNe. Since virtually 
all massive stars are born in binaries and up to 70% of them 
will interact with their companion (Sana et al. 2012), binary 
interaction may be the top contender for how massive stars 
shed their hydrogen-rich envelopes. It is possible that binary 
interactions (and other massergon; 15-loss mechanisms) can 
result in various degrees of envelope stripping and that there 
is a broad distribution of hydrogen-rich envelope masses at 
the presupernova stage of stars of any given zero-age main 
sequence (ZAMS) mass. 

Our goal here is to study the effects of substantial mass 
loss during massive star evolution on presupemova struc¬ 
ture and the resulting SN light curve. We assume that the 
mass is lost rapidly (e.g., in an unstable mass transfer event 
or through some instability), but sufficiently slowly that the 
star can re-adjust to a new equilibrium after the mass loss 
event. Rather than attempting to self-consistently simulate 
various highly uncertain mass loss mechanisms, we instead 
conduct a controlled experiment in stellar astrophysics by sys¬ 
tematically stripping material from the envelope of a fiducial 
-^ZAMS = 15 Mq star at different points of its evolution. We 
note that Bayless et al. (2015) carried out a similar study of 
the effects of mass stripping on Type II SN light curves. How¬ 
ever, they considered a 23-Mq progenitor model and stripped 
it only at the presupernova stage. 



logio(Teff/[K]) 


Figure 1. Evolutionai'y track on the Hertzsprung-Russell diagram for the un¬ 
stripped reference 15-Mq (at ZAMS) star computed with MESA. Each mark 
corresponds to a post-MS evolutionary stage at which we strip: mSGB (red 
dot) stands for “middle of the subgiant branch” defined via Tgff = 10^ K; 
hMR (yellow triangle) stands for “half maximum radius,” and MCE (cyan 
rhombus) stands for “maximum extent of convective envelope.” While these 
three stripping points are separated only by yrs , they sample an in¬ 

teresting range of structure and an order of magnitude in envelope extent. 
See Table 1 for more quantitative information on the stripping points. Each 
stripped model corresponds to the reference model up to its stripping point 
and is stripped only once. 

Figure 1 shows the evolutionary track on the Hertzsprung- 
Russel diagram for the 15-Mq reference model. Rapid mass 
loss is most likely to occur in the post-MS evolution because 
the envelope expands and becomes only weakly bound to the 
compact core. When precisely the mass loss event occurs and 
how much mass is lost will depend on the process causing 
mass loss and possibly widely varying parameters such as the 
details of binary configuration. In order to account for our 
ignorance of the details of the mass loss event, we consider 
three points (indicated by symbols in Figure 1; see also Ta¬ 
ble 1) in the post-MS evolution of our reference model that 
probe different envelope structures and span envelope radii 
from ~ 8 Oi ?0 to ~64Oi?0: 

• mSGB series: These models are stripped at T^s = 
10^ K, which marks the middle of the subgiant branch 
(mSGB). At this point, the star’s envelope has expanded 
to a radius of 79.8 Rq. Hydrogen is burning in a shell 
and the ~5 Mq helium-rich core is inert. The envelope 
is still mostly radiative with a convective layer at mass 
coordinates 5.5 — 6.5 Mq. 

• hMR series: These models are stripped when the radius 
first surpasses the half-maximum radius (hMR) of the 
reference model, R ^ 375 Rq. At this point, the enve¬ 
lope region outside a mass coordinate of m ^ 9.5 Mq 
is convective while deeper layers are radiative. There is 
hydrogen shell burning and a small core helium burning 
region. 

• MCE series: These models are stripped when the maxi¬ 
mum radial extent of the convective envelope is reached 
(MCE, at a stellar radius of i? ~ 638 Rq). The outer 
^8Mq of the star are fully convective at this point. 
There is hydrogen shell burning and a small core he¬ 
lium burning region. 
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In each series of models, we strip in units of 1 Mq, but stop 
before we reach the hydrogen shell burning region and keep 
at least IMq of the radiative layer that surrounds it. Note that 
if we considered Roche-lobe overflow in a binary system as 
the mechanism for mass loss, then only the outer convective 
layers could be lost in an unstable Roch-lobe overflow event, 
but stripping of radiative layers would not occur dynamically 
(cf. Hilditch 2001). 

The time spanned by mSGB-hMR-MCE is only of or¬ 
der 10^ yrs, which is small compared to the full lifetime of 
the unstripped reference model (^14.13 x 10® yrs). It is, 
however, very large compared to the dynamical time and 
the thermal Kelvin-Helmholtz time of the star. The latter 
is txH ~ 3/4GM^/(i?L) (Kippenhahn et al. 2013), which 
is ~1250yrs, ~425yrs, and ~125yrs at mSGB, hMR, and 
MCE, respectively. After stripping, about ~10® yrs of evolu¬ 
tion are left until core collapse. 

3.2. MESA Simulations 

We employ MESA release version 6794 (Paxton et al. 
2011, 2013) and assume solar metallicity Zq = 0.019. 
We use the Ledoux criterion for convection and follow 
Sukhbold & Woosley (2014), who set the mixing length pa¬ 
rameter ttmit = 2.0, the overshooting parameter /ov = 
0.025, and semiconvection efficiency age = 0.01, and 
do not consider thermohaline mixing. We use the wind 
mass loss prescription of Vink et al. (2000, 2001) for the 
hot MS phase and that of de lager et al. (1988) for the 
cool giant phase, both with 77 = 1.0. We limit MESA’s 
timestep by enforcing fractional changes in structure and 
thermodynamics variables of less than 10“® per timestep 
(varcontrol_target= 10 “®) and also use a customized 
timestep control that enforces a timestep that is always 
smaller than the model’s Kelvin-Helmholtz time. We use 
mesa’s standard setting for rezoning, inesh_delta_coef f 
= 1.0, andinesh_delta_coef f_for_highT = 1.5, which 
coarsens the resolution at T > 10® K and, thus, in the core 
region, where we do not currently trust MESA results (see be¬ 
low). These standard resolution setting provide a sufficiently 
resolved envelope for our SNEC explosion and light curve 
simulations. 

Eor simplicity and speed of execution, we simulate all mod¬ 
els with mesa’s default 21-isotope nuclear reaction network 
approx21 until the onset of core collapse, which is com¬ 
monly defined as the point when the infall velocity reaches 
1000 km s“®. We note that much larger (100 — 1000 isotope) 
networks are needed for an accurate treatment of late oxy¬ 
gen burning and silicon burning and of the pre-collapse neu- 
tronization in the degenerate core (e.g., Sukhbold & Woosley 
2014; Arnett, private communication). Since the treatment 
of these late burning stages has a large effect on the core re¬ 
gion out to the carbon-oxygen - helium interface (Sukhbold 
& Woosley 2014), the core structure of our MESA models is 
not reliable. Core collapse and postbounce supernova simu¬ 
lations that focus on the explosion mechanism suggest that 
the structure in the inner 1.4 — 2.5 Mq may determine if 
neutrino-driven explosions fail or succeed (O’Connor & Ott 
2011; Ugliano et al. 2012; Ertl et al. 2015). However, for our 
explosion study with SNEC, details of the core structure are 
not essential, since the outer regions and the hydrogen-rich 
envelope determine the light curve. We also artificially launch 
explosions and introduce ®®Ni by hand (see Section 2). 


3.3. Stripping Procedure 

We first evolve copies of the unstripped reference model 
to the onset of core collapse, and to the three stripping 
points: mSGB, hMR, and MCE (see Table 1 and Fig. 1). At 
these stripping points, we then restart and use MESA’s mod¬ 
ule ad justjnass to “instantaneously” remove the specified 
amount of mass. The new smaller value of the total mass 
is reached using 80 MESA “pseudo-evolution” steps (i.e., the 
structure is evolved using timesteps, but the time coordinate 
is held constant). During each step, MESA removes the largest 
amount of mass from the envelope that it can while still find¬ 
ing a hydrostatic solution to the structure equations. In most 
cases, ^75 are sufficient to reach the desired mass and the 
last ~5 have mass loss set to zero. Both during and at the 
end of the pseudo-evolution, the structure is always in global 
hydrostatic equilibrium, therefore, when the regular evolution 
resumes, no readjustment occurs. 

We strip mass in 1 Mq steps and continue the evolution of 
each model to the onset of core collapse. We refer to the un¬ 
stripped reference model simply as “unstripped,” and name 
the stripped models according to [stripping point] .[number of 
Mq strippedjM. For example, “hMR_5M” stands for a model 
that had 5 Mq stripped at hMR. 

3.4. Resulting Presupernova Structure 

We summarize the presupernova structure of our model 
set in Table 2 and Figure 2. The unstripped reference 
model reaches the presupernova stage as a 12.28 Mq RSG 
with a hydrogen-rich envelope mass of Mh ~ 7.18Mq. 
The stripped models have envelopes of systematically lower 
mass, approximately proportional to the amount of mass re¬ 
moved. While Mh varies by more than a factor of seven 
(Mh ~ 0.31 — 0.38 for the most stripped models) within 
a model series, the final stellar radius varies only by a fac¬ 
tor of ~2. Hence, the envelopes become increasingly dilute 
(lower-density) with increasing stripping. Following the Teg 
criterion of Georgy (2012), our models with less than 6 Mq 
stripped die as RSGs while models from which we strip 6 Mq 
or 7 Mq die as yellow supergiants (YSGs). It is apparent from 
Figure 2 that the choice of stripping point has almost negligi¬ 
ble influence on the final mass and structure of the remain¬ 
ing hydrogen-rich envelope. Differences in Mh and radius 
R are generally < 5% between models from different series 
that had the same amount of mass removed. An exception are 
the most extreme JVI7 cases (7Mq of hydrogen-rich material 
stripped) that show a ~ 20 % variation in their final Mh from 
0.31Mq in model MCE_7M, 0.32Mq in model hMR_7M, 
to 0.38 Mq in model mSGB_7M. The envelope in the latter 
model temporarily becomes compact when helium ignites in 
the core, which leads to less wind mass loss after stripping. 
The radii of all _7M models at core collapse are nearly identi¬ 
cal (-550 i?Q). 

The iron core mass (—1.5 — 1.6 Mq) and density profile 
is nearly identical in all models. They reach core collapse 
at central densities in the range 0.93 — 1.48 x 10®®gcm“®. 
Since the electrons in the iron core are relativistically degen¬ 
erate and the core specific entropy and electron fraction are 
roughly the same in all models, the iron core structure is very 
similar throughout the model series. More interesting are the 
large variations in the density profiles in the silicon and oxy¬ 
gen/carbon layers above the iron core, between —1.5 Mq and 
—3.3Mq, as shown in Figure 2. The presupernova structure 
in these layers appears to be very sensitive to both the amount 
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Table 1 

Summary of stripping points. We give the the criterion defining the stripping point, the unstripped reference model age, radius (R), hydrogen-rich envelope 
mass (Mu), helium-core mass (Mne), and the maximum mass stripped (max{ AM}) for each stripping point. We generate stripped models at each stripping 
point at the timestep at which the reference model exactly meets or exceeds the stripping criterion. In the criterion for MCE, Xc is the abundance of hydrogen in 
the central computational cell; is the unweighted average convective velocity in the outermost 150 computational cells, is the unweighted average 

convective velocity in the 150 computational cells above the outermost lower boundary of a convective region. If these two differ by less than 10%, the envelope 
has roughly reached its maximum extent. 


Series Name 

Stripping Criterion 

Age [Myr] 

R[R@] 

1 

1 

A^He Wq\ 

max{AM} 

Middle SGB (mSGB) 

0 

II 

13.0263 

79.8 

10.67 

3.81 

7 

Half Maximum Radius (hMR) 

R > 375 Rq 

13.0305 

381.6 

10.62 

3.87 

7 

Maximum Extent 

of the Convective Envelope (MCE) 

Xc = 0 and 

(l,surf _ ^env )/^surf <01 
V^^conv ^conv// ‘^conv 

13.0310 

638.1 

10.61 

3.88 

7 


Table 2 

Summary of the presupernova structure of the MESA models. Mp^e—SN is the total presupernova mass, Mh is the mass of the hydrogen-rich envelope, and 
Mjje, Mqo, and Mp^ are the helium, carbon/oxygen, and iron core masses, respectively. We place the iron core boundary at the first location going inward 
where the electron fraction Ye < 0.49, following the definition of Dessart et al. (2013). Eb is the binding energy of the material outside 1.4 Mq (the mass 
coordinate of the inner boundary in our SNEC explosion models) given in units of Bethe, 1B = 10®^ erg. $ 2 ^ 5 ” ^ is the compactness parameter of O’Connor 
& Ott (2011) and R, L, and Tgfr are the presupemova stellar radius, luminosity, and effective temperature, tgg is the time of shock breakout in hours after the 
onset of the thermal bomb. 


Model 

Afpre-SN [A^©] 

Mh [Mq] 

AThb [AT©] 

0 ^ 

0 

1 

ATpe [AT©] 

\Eh\ [B] 

j^pre —biN 

S2.5 

R[Rq] 

L[Lq] 

CO 

0 

tSB [h] 

unstripped 

12.28 

7.18 

5.10 

3.27 

1.51 

0.641 

0.103 

1039 

120309 

3.337 

48.52 

mSGBMM 

11.27 

6.18 

5.09 

3.28 

1.49 

0.697 

0.125 

1031 

121084 

3.355 

45.71 

mSGB_2M 

10.25 

5.16 

5.09 

3.26 

1.49 

0.617 

0.142 

1013 

119370 

3.373 

42.79 

mSGB_3M 

9.17 

4.06 

5.11 

3.27 

1.49 

0.586 

0.127 

991 

121536 

3.425 

38.74 

mSGB_4M 

7.87 

2.67 

5.20 

3.32 

1.58 

0.749 

0.138 

932 

122270 

3.537 

31.99 

mSGB_5M 

6.82 

1.61 

5.21 

3.33 

1.56 

0.734 

0.171 

828 

122984 

3.759 

25.11 

mSGB_6M 

5.94 

0.74 

5.20 

3.32 

1.54 

0.650 

0.114 

663 

123258 

4.204 

16.75 

mSGB_7M 

5.59 

0.38 

5.21 

3.33 

1.50 

0.625 

0.089 

555 

118763 

4.553 

12.25 

hMR_lM 

11.27 

6.18 

5.09 

3.28 

1.49 

0.697 

0.125 

1031 

121084 

3.355 

45.71 

hMR_2M 

10.25 

5.16 

5.09 

3.26 

1.49 

0.617 

0.142 

1013 

119370 

3.373 

42.79 

hMR_3M 

9.17 

4.06 

5.11 

3.27 

1.49 

0.586 

0.127 

991 

121536 

3.425 

38.74 

hMR_4M 

7.87 

2.67 

5.20 

3.32 

1.58 

0.749 

0.138 

932 

122270 

3.537 

31.99 

hMR_5M 

6.87 

1.68 

5.19 

3.32 

1.53 

0.658 

0.118 

843 

122179 

3.719 

25.58 

hMR_6M 

5.96 

0.77 

5.18 

3.31 

1.58 

0.709 

0.122 

676 

122065 

4.153 

17.19 

hMR_7M 

5.52 

0.32 

5.21 

3.30 

1.60 

0.706 

0.110 

551 

122645 

4.604 

11.83 

MCEMM 

11.27 

6.17 

5.10 

3.27 

1.58 

0.765 

0.102 

1032 

118857 

3.339 

45.62 

MCE_2M 

10.25 

5.16 

5.09 

3.27 

1.54 

0.644 

0.134 

1016 

120982 

3.379 

42.79 

MCE_3M 

9.17 

4.06 

5.11 

3.27 

1.53 

0.696 

0.159 

989 

119197 

3.413 

38.65 

MCE_4M 

7.88 

2.69 

5.19 

3.32 

1.58 

0.592 

0.130 

932 

122808 

3.541 

32.36 

MCE_5M 

6.87 

1.68 

5.19 

3.32 

1.52 

0.708 

0.131 

843 

121709 

3.715 

25.55 

MCE_6M 

5.96 

0.78 

5.18 

3.31 

1.55 

0.694 

0.153 

675 

122791 

4.162 

17.20 

MCE_7M 

5.52 

0.31 

5.21 

3.31 

1.56 

0.718 

0.123 

552 

122631 

4.602 

11.80 


of envelope mass stripped and the stripping point in the evo¬ 
lution. However, there are no identifiable trends that could 
be linked to amount of mass stripped and stripping point. 
Sukhbold & Woosley (2014) pointed out that the structure in 
these layers is sensitive to the treatment of (1) nuclear reac¬ 
tions and weak interactions (neutrino cooling, neutronization) 
and (2) mixing and overshooting. Both (1) and (2), in turn, 
influence the number and extent of convective shell burning 
episodes/regions in the silicon and carbon/oxygen layers. Our 
results indicate that variations in the time and amount of mass 
loss can also influence this part of presupernova stellar struc¬ 
ture. The density distribution in the affected regions deter¬ 
mines the compactness parameter of O’Connor & Ott (201 1), 

c M/Mq 

“ i?(M)/1000km ’ ^ ’ 

with the commonly adopted reference value M = 2.5Mq. 
Multiple studies (e.g., Ertl et al. 2015; Ugliano et al. 2012; 


O’Connor & Ott 2011) have demonstrated that the compact¬ 
ness parameter ^ 2.5 is a useful quantity to (at a roughly “first 
order” level) judge whether a given star is more likely to ex¬ 
plode in a supernova or collapse to a black hole without ex¬ 
plosion. Hence, the dependence of ^ 2.5 on mass loss (both 
rapid and due to winds) deserves further investigation in fu¬ 
ture work with a version of MESA with a much larger nuclear 
reaction network and a more reliable treatment of the final 
stages of stellar evolution. 

4. LIGHT CURVES 

We next present our study of the light curves from explod¬ 
ing the stripped MESA models described above. We begin by 
summarizing our basic setup in Section 4.1 . We then present 
the light curve of an unstripped model in Section 4.2 that will 
serve as a reference for the subsequent discussion of stripped 
models. We also include a discussion of how details such as 
the mixing, the nickel distribution, and how the explosion is 
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Figure 2. Density profiles as a function of enclosed mass at the onset of 
core collapse of the entire set of MESA presupernova models computed in 
this study. The solid curves include the unstripped reference model and the 
models stripped at mSGB, dashed curves show models stripped at hMR, and 
dot-dashed curves represent models stripped at MCE (cf. Figure 1 and the 
text for a discussion of these stripping points). More stripped models have 
more tenuous envelopes, but the time of stripping has negligible influence on 
the envelope structure. However, both time of stripping and the amount of 
mass stripped influence the structure of the silicon and carbon/oxygen lay¬ 
ers ai'ound 1.5 — 3.3 Mq in a complex and not obviously systematic way. 
The structure in this region may determine the outcome of core collapse 
(cf. O’Connor & Ott 2011; Ugliano et al. 2012; Ertl et al. 2015). 


initiated impact the light curves in Section 4.3. Finally, we 
present our main study of the light curves from models with 
varying levels of stripping in Section 4.4. 

4.1. Explosion and Light Curve Setup 

In all explosion models, we excise the inner 1.4 Mq, as¬ 
suming that this part collapsed and formed a neutron star. 
We then map (via linear interpolation) the hydrodynamic and 
compositional variables from MESA to a 1000 cell grid in 
SNEC. We choose the grid spacing so that resolution is con¬ 
centrated in the interior, where the thermal bomb is placed, 
and near the surface, where the photosphere is initially lo¬ 
cated. In our fiducial resolution calculations, the innermost 
cell has a mass of Am = 6.5 x 10“® Mq and the surface cell 
has a mass of Am = 6.5 x 10“® Mq. The lowest resolu¬ 
tion in our fiducial setup is Am = 0.065 Mq at mass coordi¬ 



Figure 3. Mass fractions of the key elements hydrogen, helium, oxygen, 
and carbon in the unstripped reference model before (upper panel) and af¬ 
ter (lower panel) boxcar smoothing. Note that the smoothed profile here is 
shown for a grid with uniform spacing in mass. In our production setup that 
uses a non-uniform grid with higher resolution in the innermost and outer¬ 
most regions, small jumps in the composition profiles remain, but have no 
influence on the light curve. We excise the mass inside the shaded regions 
before launching the explosion. 


nates between ~2.5 and ~5 Mq (at around grid cell 100). The 
mass of cells between the innermost cell and the coarsest cell 
changes according to geometric progression with a fixed ratio 
between two consecutive cells > 1. Between the coarsest cell 
and the surface cells we refine by geometric progression with 
a fixed ratio between two consecutive cells < 1. Examples of 
the grids may be found in Appendix A. The release version of 
SNEC contains a routine to generate a variety of grid setups. 

The explosion is initiated by applying a “thermal bomb” to 
the innermost region of the model, just above the mass cut 
in a way similar to previous work (e.g., Blinnikov &. Bar- 
tunov 1993; Bersten et al. 2011; Aufderheide et al. 1991). 
The energy of the bomb is added to the right-hand-side of 
Equation (2) during a chosen time interval and in a range of 
mass from the inner boundary, both exponentially attenuated. 
Eor our fiducial model, summarized in Section 4.2, the bomb 
is spread over 0.02 Mq and injected in Isec. The thermal 
bomb mechanism, implemented this way, typically gives a 
few percent more energy to the system than dialed-in, due to 
the smooth exponential attenuation. This small excess in en¬ 
ergy is recorded and accounted for in SNEC’s global energy 
balance. We find that SNEC conserves total energy to better 
than 1% in a full-physics model followed to 150 days past 
explosion. 

An alternative way of phenomenologically modeling an SN 
explosion is the “piston mechanism,” which has been used, 
for example, in the work by Eastman et al. (1994), Utrobin 
(2007), Kasen & Woosley (2009) and Dessart et al. (2013). 
Eor the same amount of injected energy, the thermal bomb 
and piston mechanism give nearly identical light curves, as 
was discussed in Bersten et al. (2011) and confirmed by our 
own simulations. However, when a reaction network is in¬ 
cluded, piston and thermal bomb may result in different nu- 
cleosynthetic yields (Young &. Eryer 2007). We have imple¬ 
mented a piston in SNEC, but do not use it in this work, since 
the thermal bomb makes it easier to control the energy of the 
explosion. 

It has long been realized in SN light curve modeling that 
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Figure 4. Bolometric light curve of the unstripped reference model. Time 
is given relative to the onset of the thermal bomb driving the explosion. 
Shock breakout occurs at day 2.03 through the reference model’s surface 
at ~1039 Rq. The black graph is the fiducial light curve obtained with our 
standard parameter choices, including boxcar smoothing as an approxima¬ 
tion of mixing during the explosion (cf. Figure 3 and §4.1). The red graph 
represents the unmixed case with steep compositional gradients (top panel of 
Figure 3). The inset plot shows shock breakout and the very early light curve. 
We note that during shock breakout the photosphere is located in the outer¬ 
most cell of SNEC’s grid and spatially poorly resolved. Thus the light curve 
predicted by SNEC around the time of shock breakout is likely not reliable 
(cf. Ensman & Burrows 1992). 

sharp gradients in the composition prohles of the progenitors 
may result in artihcial light curve features that are not ob¬ 
served in real SNe. For example, Utrobin (2007) points out 
a pronounced bump at the end of the plateau, followed by an 
abrupt decrease of the bolometric luminosity for a model with 
unmixed chemical composition (see their Figure 16). Two- 
and three-dimensional simulations of shock propagation in 
core-collapse SN explosions (e.g., Kifonidis et al. 2003,2006; 
Wongwathanarat et al. 2015) show that effective mixing oc¬ 
curs during the shock propagation through the progenitor due 
to the Rayleigh-Taylor and Richtmyer-Meshkov instabilities. 
In this process, hydrogen and helium get mixed into the inner 
layers, while metal-rich clumps, and, in particular, ®®Ni, may 
penetrate outwards up to ^3000 km and more in the ve¬ 
locity prohle. Lacking a physical mechanism for the mixing 
in our one-dimensional code, we apply an artihcial “boxcar” 
averaging, as used, for example, in Kasen & Woosley (2009) 
and Dessart et al. (2012, 2013). We run a boxcar with a width 
of 0.4 Mq through the model four times until we obtain a 
smooth prohle (details of this procedure are available in the 
SNEC notes document on the SNEC website). As an example. 
Figure 3 depicts the non-mixed (top panel) and mixed (bot¬ 
tom panel) mass fractions of hydrogen, helium, oxygen, and 
carbon in the unstripped reference model. 

Finally, we assume a hxed amount of ®®Ni of 0.05 Mq in all 
our models, which is roughly the average amount deduced for 
SNe IIP (which have a range of ~ 0.01—0.1 Mq, see Kasen & 
Woosley 2009; Smartt et al. 2009). We distribute it uniformly 
in the interval between the excised mass of 1.4 Mq and some 
chosen mass coordinate (5 Mq for the reference runs, which 
is near the edge of the helium core) at the expense of other 
elements before smoothening the composition. 

4.2. Fiducial Light Curve of the Unstripped Reference 
Progenitor Model 

Figure 4 shows the light curve of the unstripped reference 
model. The explosion was initiated at f = 0 with a ther- 
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Figure 5. Bolometric light curves of the unstripped reference progenitor 
computed with different mass range over which the thermal bomb is spread 
(top panel) and different durations of the thermal bomb (bottom panel). All 
other parameters are those laid out in §4.1. Time is relative to the onset of the 
thermal bomb. 
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Figure 6. Bolometric light curves for the unstripped reference progenitor 
model with different degree and extent of initial ®®Ni mixing and all other 
parameters as laid out in §4. 1. The helium core mass is 5.1 Mq and if Ni is 
mixed smoothly into the hydrogen-rich envelope (green graph), then the light 
curve’s “knee” feature visible when the photosphere drops into the helium 
core disappears. 

mal bomb resulting in an asymptotic kinetic energy of E = 
lO^i gj.g (as described in §4.1). We added Mn; = 0.05 Mq, 
mixed up to a mass coordinate of 5 Mq (but not into the hy¬ 
drogen envelope, though boxcar smoothing introduces some 
nickel into the hydrogen-rich region.). These are the fiducial 
explosion parameter choices used throughout this study. 

The light curve of the unstripped reference model shows all 
the traditional hallmarks of a typical SNe IIP (e.g., Falk & 
Arnett 1977; Eastman et al. 1994; Filippenko 1997). Shock 
breakout occurs when the optical depth is less than ~ c/u at 
a time of 2.03 days after the onset of explosion. The bolo¬ 
metric luminosity peaks at L = 3.4 x lO^^ergs”^ with an 
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effective temperature of Teff = 1-7 x 10® K. The subsequent 
cooling phase (discussed in detail in Nakar & Sari 2010) lasts 
for ~ 19 days. At this point, the ejecta have expanded and 
cooled so much that hydrogen recombination sets in (starting 
already at T ~ 7500 K) and powers the plateau phase with 
very slowly decreasing effective temperature that varies from 
^ 6000 K at ^ 35 days down to ~ 5000 K at ~ 90 days. 
The recombination wave and, consequently, the photosphere 
moves inward in mass coordinate, but due to the overall ex¬ 
pansion stays at roughly constant radius, resulting in a rela¬ 
tively small variation in luminosity during the recombination 
phase of the plateau from day ^ 19 to day ~ 80 — 90 (this 
phase is investigated analytically in Goldfriend et al. 2014). 
The slow decline that is apparent in the plateau phase shown 
in Figure 4 occurs because of the combined effects of the pho¬ 
tosphere receding slightly in radius and the effective temper¬ 
ature slowly decreasing (e.g., Eastman et al. 1994; Woosley 
1988). 

The plateau ends when the photosphere reaches the helium 
core. Helium recombines at T > 10"^ K whereas the photo- 
spheric temperature is T ^ 5000 K, recombination acceler¬ 
ates dramatically, and both the radius of the photosphere and 
the luminosity decrease rapidly. Note that it is common in the 
SN IIP theoretical light curve literature to define the plateau 
duration as the time from shock breakout to the drop when 
the photosphere reaches the helium core (Kasen & Woosley 
2009; Popov 1993). We adopt this definition of plateau dura¬ 
tion in this paper. The small “knee” or “bump” feature in the 
drop of the fiducial light curve around day 90 in Figure 4 is 
due to the additional luminosity input from radioactive ®®Co 
(from the ®®Ni—/-^^Co—>^®®Fe decay chain) that is uncovered 
as the photosphere sweeps through the helium core through¬ 
out which ®®Ni was mixed initially. This feature is sensitive 
to the degree and implementation of mixing and is unlikely to 
be robust (see next Section 4.3). Finally, the tail of the light 
curve, after day ^100, is powered exclusively by the radioac¬ 
tive decay of ®®Co. 

4.3. Sensitivity of the Fiducial Light Curve to Mixing, 
Thermal Bomb Parameters, and Nickel Distribution 

The red curve in Figure 4 highlights the effect of steep com¬ 
positional gradients on the light curve in comparison with the 
result obtained with compositional smoothing (black graph; 
“boxcar averaging”; Section 4.1) that we use to mimic multi¬ 
dimensional mixing during the explosion. If exploded with¬ 
out smoothing, hydrogen-rich material transitions discontin- 
uously to helium-rich material (cf. Figure 3), which leads to 
more rapid recombination, a more abrupt drop of the photo¬ 
sphere radius, and a steeper decline of the luminosity. Al¬ 
though observations of most SNe do not generally reveal such 
abrupt drops, some subclasses of SNe may have rapidly drop¬ 
ping photospheric velocities as discussed by Piro & Morozova 
(2014) because of this same effect of rapid helium recombi¬ 
nation. For all other light curves presented in this paper, we 
use the smoothed composition profiles. 

In Figure 5, we explore the sensitivity of the fiducial light 
curve to (top panel) variations in the amount of mass over 
which the thermal bomb is spread and (bottom panel) varia¬ 
tions in the duration over which the energy is injected. While 
the details of the energy injection will depend on the ac¬ 
tual physical explosion mechanism (e.g., Bethe 1990; Janka 
2012 ), it is reassuring that the light curve is fairly insensitive 
to both mass spread and duration of energy injection. We re¬ 
mind the reader that for the light curves shown in Sections 



Figure 7. Upper panel: Bolometric light curves from the mSGB grid of pro¬ 
genitor models. Lower panel: Light curves from the models ‘unstripped’ and 
‘mSGB_7M’. Dashed lines show the heating rate from the radioactive de¬ 
cay of ^®Ni deposited in each model after taking into account leakage of the 
gamma rays. In the most stripped model, gamma rays leak out faster than in 
the unstripped model, impacting the late time light curve. 

4.2 and 4.4, the thermal bomb is spread over 0.02 Mq and 
its duration is 1 s (in self-consistent multi-dimensional core¬ 
collapse SN simulations most of the energy injection appears 
to occur within ~ 1 — 2 s; e.g., Bruenn et al. 2013). 

Finally, Figure 6 shows the dependence on ®®Ni mixing. 
The overall effect of ®®Ni mixing is modest. The general 
trend is that with the initial mixing of ®®Ni to increasing 
mass coordinates, the contribution to the luminosity of the 
®®Ni—)-®®Co—/-^^Fe decay chain becomes more prominent. At 
later times, the luminosity due to radioactive decay becomes 
smaller with increasing mixing due to the fact that more 
gamma rays escape from the model without being absorbed 
(as also described by Young 2004; Utrobin 2007; Bersten 
et al. 201 1). Note that if ®®Ni is mixed far into the hydrogen- 
rich envelope, the “knee” feature at the end of the plateau 
disappears. This is consistent with the findings of Kasen & 
Woosley (2009), who mixed ®®Ni into the hydrogen-rich en¬ 
velope in their models. 

4.4. Light Curves as a Function of Mass Stripping 

Figure 2 demonstrates that the structure of the hydrogen- 
rich envelope and of the outer helium core are essentially in¬ 
dependent of the point at which we remove mass for the set 
of stripping points we choose in this study (cf. §3.4). This 
suggests that the resulting light curves should be independent 
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Figure 8 . The same bolometric light curves as in Figure 7, but instead fo¬ 
cusing on early times around shock breakout. Time is given relative to the 
onset of the thermal bomb driving the explosion. Solid parts of the curves 
start from the time at which the photosphere moves inwards in the grid space, 
while dashed parts indicate that the photosphere is located in the outermost 
grid cell and is spatially poorly resolved. Note that more stripped models have 
a higher luminosity in the post-breakout cooling phase and a faster evolving 
light curve. 
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Figure 9. Photospheric velocities Uph in the mSGB model series, including 
the unstripped reference model. These velocities are reliable only until the 
photospheric temperature drops below T ~ 10^'^® K below which we can¬ 
not accurately estimate the location of the photosphere due to low-T opacity 
limitations in SNEC. We indicate the unreliable part of Upjj by plotting it 
in dashed lines. Note that the temperature drops most rapidly in the most 
stripped models that lack a plateau phase (cf. Figure 7). 


of the stripping point and we check this assertion later in this 
section. Here, we operate under the assumption that it is true 
and focus our discussion on the model series stripped at the 
middle of the suhgiant branch (mSGB). 

In the top panel of Figure 7, we show bolometric mSGB 
series light curves obtained for a final ejecta kinetic energy 
of E'kin = 10®^ erg and Mni = 0.05 Mq (and all other ex¬ 
plosion parameters as specified in §4.1). The early, shock 
breakout and cooling part of the light curves is shown in Fig¬ 
ure 8. Shock breakout itself is not well resolved (i.e., the 
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Figure 10. Bolometric light curves for all models from the three different 
MESA series, described in the Tables 1 and 2. The explosion setup is identical 
in all cases and as described in §4.1. All models are stripped after they have 
left the main sequence, but the precise point of stripping has little influence 
on the resulting light curve. 


photosphere is in the outermost grid cell) and thus the light 
curves in this phase are unreliable (cf. Ensman & Burrows 
1992). Once the photosphere begins to move inward into the 
expanding envelope, the light curves predicted by SNEC be¬ 
come robust. Models with greater amounts of mass stripped 
have higher luminosities in the cooling phase and decay more 
rapidly. In more stripped models, the SN shock has to prop¬ 
agate through less envelope mass that is more tenuously dis¬ 
tributed. This results in a higher shock velocity at early times 
as shown in Figure 9 (and earlier breakout; as shown in Fig¬ 
ure 8), which leads to a hotter photosphere and a more rapid 
expansion of the ejecta. This translates directly to a higher 
initial luminosity and a more rapid decay of the light curve. 

From the top panel of Figure 7 we see that as long as there 
is a substantial amount of hydrogen-rich material left, a clear 
(if very short) plateau due to recombination can be made out. 
In our mSGB model series, this is until models mSGB_4M 
and mSGB_5M, from which we strip AMq and 5 Mq and 
which have 2.67 Mq and 1.61 Mq of hydrogen-rich material 
left, respectively. The two most stripped models of this se¬ 
ries, models mSGB_6M and mSGB_7M, have only 0.74 Mq 
and 0.38 Mq of very tenuous hydrogen-rich envelope left, re¬ 
spectively. This does not appear to be sufficient to lead to 
any plateau and the photosphere recedes very quickly in these 
models. Instead, models mSGB_6M and mSGB_7M show 
a clear peak between ~ 20 — 30 days that is analogous to 
the nickel-powered peak that is seen in all types of hydrogen- 
dehcient (i.e.. Type I) SNe. 

Figure 7 shows that stripping of hydrogen-rich envelope 
mass also has an effect on the late-time radioactively pow¬ 
ered part of the light curve. The late-time light curves exhibit 
changes with mass stripping because there is earlier leakage 
of gamma rays from the more highly stripped models. This is 
shown by the lower panel of Figure 7. It shows the unstripped 
reference model and the most stripped model (mSGB_7MQ) 
in comparison with the total amount of heating deposited in 
each model due to the radioactive ®®Ni—>®®Co—>®®Fe decay 
chain. The most stripped model has less radioactive heating 
at late times. 
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Figure 10 compares the bolometric light curves of mod¬ 
els stripped at the middle of the subgiant branch (mSGB, the 
model series we focus on), half-maximum radius (hMR), and 
maximum radial extent of the convective envelope (MCE); 
cf. §3.1 and Table 1. The light curves of models with the same 
amount of mass stripped at different times are nearly identical, 
supporting our initial assertion on the basis of Figure 2. Any 
variation between light curves of models stripped at mSGB, 
hMR, and MCE is arguably smaller than the level of system¬ 
atic uncertainty inherent to SNEC’s approximate (i.e., equi¬ 
librium diffusion) way of predicting these light curves. The 
robustness of the light curves suggests that comparisons with 
observations may allow reliable conclusions about the amount 
of hydrogen left after (rapid) mass loss events. However, we 
note that we exclusively consider post-MS mass loss. Early 
large-scale mass loss events on the MS may lead to different 
outcomes, which should be explored in future work. 

In Figure 11, we plot approximate absolute magnitudes 
of the mSGB series light curves in the IRVB- and t/-bands 
and in Figure 12, we focus on the V-band. We obtain the 
band light curves by assuming black body emission from the 
photosphere and using the bolometric corrections from Ofek 
(2014). When interpreting these light curves, one should keep 
in mind two important caveats: ( 1 ) When the whole ejecta 
becomes optically thin, the luminosity has a large contribu¬ 
tion of ®®Ni/^®Co from above the photosphere. For this rea¬ 
son, we terminate the curves at the points where the lumi¬ 
nosity contribution due to ^®Ni/®®Co above the photosphere 
amounts to more than 5% of the total luminosity. (2) As 
was demonstrated in Kasen & Woosley (2009), the U- and 
Z?-bands of the light curves cannot be adequately reproduced 
by a one-temperature equilibrium-diffusion code like SNEC, 
because these bands are strongly influenced by iron group 
line blanketing after a few tens of days (see, e.g.. Figure 8 
of Kasen & Woosley 2009). This causes a much faster de¬ 
cline of the U- and B-band light curves. However, the IR- and 
V-bands are still similar to a single temperature black body 
spectrum, and thus these bands are more accurately captured 
by SNEC. 

Finally, we study the sensitivity of our light curves to dou¬ 
bling the explosion energy to 2 x 10^^ erg and to doubling the 
amount of initially present ^®Ni to 0.1 Mq in Figure 13. The 
qualitative light curve changes are overall as expected from 
previous work (e.g.. Young 2004; Utrobin 2007; Kasen & 
Woosley 2009): More energetic explosions have brighter, but 
faster evolving light curves and an increased amount of ^®Ni 
prolongs the plateau and results in a higher late-time luminos¬ 
ity. Increasing the amount of ®®Ni also results in a more pro¬ 
nounced second light curve peak in the most-stripped mod¬ 
els mSGB_6M and mSGB_7M. We summarize the connec¬ 
tion between these calculations and observed SN light curves 
in the discussion below. 

5. DISCUSSION AND CONCFUSIONS 

We presented the new open-source SuperNova Explosion 
Code (SNEC) for investigating supernova (SN) explosions 
and the resulting light curves. SNEC, while currently limited 
to equilibrium (single temperature) radiative diffusion, is the 
first such code that is publicly available and this paper will 
serve as a reference and starting point as SNEC is utilized and 
improved by both us and the broader community in the future. 

As a hrst application of SNEC, we studied the explosions 
of Mzams = 15 Mq stars with varying levels of rapid post- 
MS mass loss at different times in the stars’ evolution from 


the main sequence to the supergiant stage. We evolved these 
stars to the onset of core collapse with the open-source MESA 
stellar evolution code and then exploded them with SNEC us¬ 
ing a thermal bomb resulting in an asymptotic explosion en¬ 
ergy of 10^^ erg. At three different times during the evolution 
to the supergiant stage, we systematically stripped hydrogen- 
rich material in units of 1 Mq, leaving in the most extreme 
case only a thin radiative hydrogen-rich layer above the hy¬ 
drogen shell burning zone. In this experiment in massive star 
evolution, we find that the time of stripping has essentially no 
influence on the structure of the envelope and thus most of the 
SN light curve. Stars with more than 1.5 — 2 Mq of hydrogen- 
rich material left die as red supergiants with R > 900 Rq and 
our most stripped star (0.31 Mq of hydrogen-rich envelope 
left) dies as a yellow supergiant with a still extended, very 
tenuous envelope of i? ~ 550 Rq. 

We find that the time of stripping and the amount of mass 
stripped has a big but not clearly systematic effect on the 
structure of the layers immediately surrounding the iron core 
(mass coordinate ~1.5 —3.3 Mq). The structure in this partic¬ 
ular region has been shown to be highly relevant for deciding 
the ultimate outcome of core collapse (explosion/no explo¬ 
sion, black hole/neutron star remnant; O’Connor & Ott 2011; 
Ugliano et al. 2012; Ertl et al. 2015). Our results suggest the 
need for a detailed study of the sensitivity of presupernova 
stellar structure to large amounts of rapid mass loss (e.g., via 
unstable mass transfer in a binary). 

The light curves resulting from our set of presupernova 
models show SN IlP-like morphology for models with more 
than ^1.5 — 2 Mq of hydrogen-rich envelope material left at 
the presupernova stage. The most stripped models (<1.5 Mq 
of hydrogen-rich envelope left) have higher luminosity in the 
post-breakout cooling phase, but show no plateau, but a sec¬ 
ond peak around 20 — 30 days due to energy input from the ra¬ 
dioactive decay of ®®Ni/®®Co that is uncovered by the rapidly 
receding photosphere in these models. 

In those models that show a plateau in their bolometric light 
curves, the duration of the plateau phase varies in the range 
^ 20 — 100 days (we include both the cooling and the recom¬ 
bination phases of the plateau in the plateau duration; Kasen 
& Woosley 2009; Popov 1993), with plateau length decreas¬ 
ing with decreasing mass of hydrogen-rich envelope material. 
In nature, most SNe IIP show plateaus of ^ 80 — 100 days 
(e.g., Poznanski et al. 2009; Arcavi et al. 2012), although there 
may be some evidence for a subset of shorter plateaus (Ander¬ 
son et al. 2014). The completely mass stripped SNe Ib/c show 
no plateau at all. Both SNe IIP and SNe Ib/c are relevant to 
our study given the inference that many SN IIP progenitors 
have ZAMS masses around ~ 15 Mq (Smartt et al. 2009) 
and arguments that most SNe Ib/c must come from a similar 
mass range (Smith et al. 2011). The apparent paucity of ob¬ 
served short and intermediate-length plateaus suggests that in 
nature hydrogen mass loss is an all or nothing process, at least 
for the ZAMS mass we consider. This is perhaps not surpris¬ 
ing given that for Mzams < 20 Mq radiative driven winds 
are rather weak in normal prescriptions and appreciable mass 
loss can probably only occur from events like binary interac¬ 
tions, which would not be expected to, say for example, rip 
off just ^ 50% of the mass. 

Our most stripped models still have « 0.3 — 0.4 Mq of hy¬ 
drogen present and thus should make some connection with 
SNe Ilb. Indeed, the light curves of these progenitors show 
two distinct peaks, similar to the morphology of many SN lib, 
where the hrst peak comes from the shock cooling of the re- 
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Figure 11. Light cui'ves in absolute magnitude in IRVB and U bands obtained with SNEC for our mSGB model set. The time is given relative to the onset of 
energy injection by the thermal bomb. The light curves start when the photosphere no longer coincides with the outermost grid cell in the SNEC calculation. 
Shock breakout, which is in the UV, would be visible in U band, but is not shown. The curves are terminated at the point at which the explosion begins to 
transition to the nebular phase and the black body approximation underlying the band light curves is no longer valid (we define this point at the time at which 5% 
of the luminosity comes from above the photosphere due to gamma-ray deposition). We note that real SN light curves fade away in U and B bands considerably 
faster than predicted by SNEC (e.g., Kasen & Woosley 2009; Dessart et al. 2013). This is due to line blanketing by iron group elements that is unaccounted for 
in our models. 



Figure 12. Light curves in absolute V-band magnitude obtained with SNEC 
for the mSGB model set. The time is given relative to the onset of energy 
injection by the thermal bomb. The V-band light curves of more stripped 
models evolve (rise, decay) faster than those of less stripped models. 


maining surface hydrogen and the second from radioactive 
heating (Woosley et al. 1994; Bersten et al. 2012; Nakar & 
Piro 2014). However, in detail, the width of the first peak is 
too large, which is likely due to our models having too much 
hydrogen still present (Nakar & Piro 2014). In addition, the 
second peak in our light curves is sometimes too dim in com¬ 


parison to observed SNe Ilb due to our models having some¬ 
what less ®®Ni (Lyman et al. 2014). However, there is a lot 
of diversity and uncertainty in the amount of ^®Ni produced 
in SNe lib (see the work and discussions in Shigeyama et al. 
1994; Bersten et al. 2012; Ergon et al. 2015), thus we plan to 
investigate this in more detail in future work. 

Another interesting connection to consider is how our re¬ 
sults relate to SNe IIL. There has been a long-standing discus¬ 
sion on whether they form a continuous sequence of events 
that smoothly transition to SNe IIP. The idea that SNe IIL 
are instead distinct from SNe IIP has been argued for by Ar- 
cavi et al. (2012) and Faran et al. (2014), but more recently it 
has been shown that SNe IIL actually show a significant drop 
in their light curves at late times (^100 days, Valenti et al. 
2015), much like SNe IIP. If there was a continuous range 
of events from SNe IIL to SNe IIP, then a natural physical 
mechanism to consider is gradual loss of the outer hydrogen, 
where SNe IIL would be on the hydrogen-poor side. At least 
for the ZAMS mass we consider here (ISM©), this does not 
appear to be the case. First, our results show that intermediate 
levels of hydrogen mass loss simply shorten the plateau length 
which is different from SNe IIL, which appear to have roughly 
normal duration, but steeply declining “plateaus” when they 
are followed for a sufficient amount of time (Valenti et al. 
2015). Second, SNe IIL are on average more luminous than 
SNe IIP by ^ 1.5 mag in the optical during the first ~ 10 days 
(Patat et al. 1993, 1994; Anderson et al. 2014; Faran et al. 
2014; Sanders et al. 2015). Our more stripped models show 
slightly higher luminosities at early times, but not nearly ex- 
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Figure 13. Bolometric light curves for the mSGB model set as in Figure 7, 
but with variations in explosion energy and ®®Ni mass. The time is given rela¬ 
tive to the onset of energy injection by the thermal bomb. Top panel: Compar¬ 
ison of light curves of models with the hducial hnal kinetic energy (10®^ erg, 
dashed curves) and with twice that energy (2 X 10®^ erg, solid curves). In¬ 
creasing the explosion energy leads to brighter, more rapidly evolving explo¬ 
sions in agreement with previous work. Bottom panel: Comparison of light 
curves of models with the hducial ®®Ni mass (0.05 Mq, dashed curves) and 
with twice that amount of ®®Ni (0.1 Mq, solid curves). More nickel leads to 
extended plateaus and brighter radioactive tails. The qualitative changes due 
to variations in explosion energy and ®®Ni mass are in agreement with what 
was found in previous work (e.g., Kasen & Woosley 2009; Young 2004). 

treme enough. Overall, however, our findings, in combina¬ 
tions with recent observations (Valenti et al. 2015), appear to 
argue that perhaps SNe IIL do not necessarily have less hydro¬ 
gen, but the hydrogen mass is distributed in a different way. 
The brighter early light curves would argue that SNe IIL have 
material at a larger radius (Piro & Nakar 2013). The occur¬ 
rence of narrow line features in SNe Iln that might otherwise 
look somewhat like an SN IIL (Smith et al. 2015) might argue 
for some contribution from circumstellar material. SNEC is 
well-suited for addressing these ideas in a systematic way in 
future work since various mass and density distributions can 
easily be implemented to investigate what is in fact needed to 
reproduce SN IIL light curves. 

Future work will be directed toward exploiting SNEC’s 
current capabilities for the systematic and reproducible light 
curve modeling for a broad range of SN explosions, but also 
toward improving SNEC’s transport solver and opacity mi¬ 
crophysics. In a first step, we will upgrade SNEC to handle 
separate radiation and matter temperatures with the long-term 
goal of constructing an open-source multi-group radiation- 


hydrodynamics code. Of course, input from the community 
will be especially critical for steering SNEC’s further evolu¬ 
tion and we look forward to the community’s feedback. 
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APPENDIX 

A. COMPARISON WITH BERSTEN ET AL. (2011) 

The work of Bersten et al. (2011) is among the codes that are most similar to SNEC, so it is important to compare light curves. 
This is done using their progenitor model for SN 1999em, kindly provided by the authors (details of the structure and composition 
of the model may be found in their work). We use the same explosion energy of E = 1.25 x 10®^ erg. We also use a step function 
for the opacity floor, i.e., 0.01 cm^ for material with Z < 0.3 and 0.24 cm^ g“^ for material with Z > 0.3. This is the only 
model in this paper where we use this opacity law, since we want to explicitly follow Bersten et al. (201 1). 

The left panel of the Eigure 14 shows bolometric light curves generated with SNEC and the light curve taken from Eigure 5 of 
Bersten et al. (201 1), together with the observational data for SN 1999em. The green light curve was generated with SNEC using 
the original 200-cell grid setup of the model of Bersten et al. (2011). Two other SNEC light curves, which were generated from the 
same model mapped onto grids with 1000 and 3000 cells, disagree with the 200 grid cells curve at early times as well as during 
the transition between plateau and ®®Ni tail. We have identified two things that may be key contributors to these differences. 

Eirst, the differences at early times before the plateau phase may be explained by differences in resolution and the resolution 
gradient between the 200-cell grid on the one hand, and the 1000-cell and 3000-cell grids on the other hand. To further illustrate 
the differences in resolution, we plot in Eigure 15 the mass resolution as a function of enclosed mass for the three grids. Note that 
during the first 50 days of the light curve, the photosphere moves inwards from a mass coordinate of 19 Mq (the total mass of 
the model) to a mass coordinate of 16 Mq. The 16 — 19 Mq region is precisely where Eigure 15 shows very large differences in 
resolution between the different grids. The original 200-cell grid of the double poly tropic model of Bersten et al. (2011) has the 
finest resolution near the surface of the star, but very rapidly changes to the coarsest resolution in the bulk of the model. In fact, 
we encountered numerical difficulties exploding the 200-cell model. We find that we have to use an outer boundary condition 
that is different from that described in Section 2.1 of this paper; for numerical stability, we find that we have to impose zero 
temperature at the surface of the star instead of constant luminosity in the two outermost grid points. The numerical evolutions 
and light curves from the 1000-cell and 3000-cell runs, however, have negligible dependence on the switch between these two 
boundary conditions. Note that the light curves from the 1000-cell and 3000-cell runs lie on top of each other, demonstrating that 
our SNEC results are numerically converged. 

Second, the difference in the transition from the plateau to the ®®Ni tail may be explained by the difference in opacities we 
use for the 200-cell run on one hand, and 1000-cell and 3000-cell runs on the other hand. Eor the 200-cell run, we use opacity 
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Figure 14. Left panel: A comparison of bolometric light curves for SN 1999em. The black and red lines show SNEC light curves generated with 1000 and 3000 
grid cells, respectively. The blue circles show observational data and the blue curve is the light curve of Bersten et al. (201 1), both taken from their Figure 5. The 
green graph shows the closest light curve we could obtain to the one of Bersten et al. (201 1) (see text for details). The dashed line gives the total power input due 
to radioactive decay. Right panel: Photospheric velocity, calculated with SNEC for the same model and resolutions as in the left panel. The blue circles show 
observational data and the blue graph shows the result of Bersten et al. (201 1), both taken from their Figure 6. Note, that we do not take the opacity floor into 
account when determining the location of the photosphere, as described in Section 2.2. 


tables that are different from those described in Section 2.2: In the low-temperature region (10^ < T < K) where 

OPAL tables are not available, we employ the Ferguson et al. tables (Ferguson et al. 2005) for all densities, temperatures, and 
compositions. These tables depend on hydrogen mass fraction, density, and temperature, and otherwise assume simply rescaled 
solar composition. Hence, we ignore the dependence of the opacity on variations of the carbon and oxygen mass fractions. We 
find that this approach is essential for reproducing the light curve of Bersten et al. (201 1). Note, however, that these authors do 
not explicitly state how they treat the opacity in low-temperature, Carbon/Oxygen rich regions. The SNEC light curves generated 
with 1000 and 3000 grid cells use our standard opacities as described in Section 2.2. As mentioned in that section, the choice 
of opacity in the low-temperature, high carbon and oxygen mass fractions region weakly influences the evolution of the system, 
because most of these opacities lie below the opacity floor. However, it does impact the position of the photosphere during the 
transition between the plateau and the ®®Ni tail, as can be seen in Figure 14. We emphasize that the green light curve in the left 
panel of Figure 14 cannot be reproduced with the standard version of SNEC. We provide it here only to demonstrate how closely 
we can approach the results of Bersten et al. (201 1). Nevertheless, this shows that overall we find reasonable agreement between 
the light curves shown in Figure 14 and have an understanding of the small differences. 

We also compare the velocity evolution calculated with the two codes. The right panel of Figure 14 shows the expansion 
velocity at the position of the photosphere, calculated with SNEC using the same progenitor model and the grid resolutions as in 
the left panel of Figure 14. The observational data are taken from Bersten et al. (201 1). As in the left panel of Figure 14, the 200 
grid cells curve shows the best agreement with Bersten et al. (201 1). 

As an additional point, it is important to understand whether our treatment applies to light curves at the earliest times. The 
typical grid setup in SNEC focuses resolution close to the surface of the progenitor star in order to ensure the photospheric 
region is well resolved as early as possible (as discussed in Section 4.1). A detailed analysis of the very early light curve around 
shock breakout was carried out by Ensman & Burrows (1992). In particular, these authors compared the results obtained with 
a two-temperature radiation-hydrodynamics treatment and with a one-temperature flux-limited equilibrium diffusion treatment 
similar to SNEC. They concluded that the two approaches give consistent results for light curve and photospheric temperature, 
provided the surface grid resolution is sufficiently fine so that at the onset of shock breakout a few tens of grid cells are covering 
the optically thin region outside the photosphere. In our SNEC models, we find that this level of photosphere resolution is not 
practical for the large model grid we are considering. However, even with our standard gridding, we find that after the first few 
hours of the explosion that this is sufficient to resolve the photosphere and produce a reliable light curve in the cooling phase 
after shock breakout. 

In order to demonstrate this point, in Figure 16 we plot the early light curves around shock breakout obtained with SNEC 
using 1000 and 3000 grid cells and the same grid setup as described above for the same progenitor model as in Figure 14. The 
model with 3000 grid cells always has at least a few grid cells above the photosphere, whereas in the lower-resolution model, the 
photosphere is remains in the outermost grid cell until it becomes resolved only about one day after breakout. Once this happens, 
the results of the 1000 and 3000 agree nearly perfectly. But even at earlier times, the two light curves agree surprisingly well. 
This can be understood from the fact that already shortly after shock breakout the light curve is determined by the diffusion of 
light from the shock heated envelope and should not depend anymore on the resolution of the surface layers. 
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Figure 15. Mass resolution as a function of mass coordinate for the 200, 1000 and 3000 grid cells simulations from Figure 14. The grid having 200 cells is the 
original grid of the model, which we received from the authors of Bersten et al. (201 1). 



Figure 16. Bolometric light curves around the time of shock breakout using the SN 1999em progenitor model. The black solid and blue dashed curves are SNEC 
light curves generated with 1000 and 3000 cells, respectively. The higher-resolution model always resolves the photosphere, and although the lower-resolution 
model does not resolve it until ~2.2 days, it agrees quite well. 


B. COMPARISON WITH DESSART ET AL. (2013) 

As an additional comparison, we consider the light curves obtained by Dessart et al. (2013) from a MESA presupemova pro¬ 
genitor. This is exploded with the Lagrangian ID hydrodynamics code VID (Livne 1993; Dessart et al. 2010), which is then 
mapped at day 10 to their non-LTE radiative transfer code CMFGEN (Hillier & Dessart 2012) that assumes homologous expan¬ 
sion. CMFGEN provides a more detailed treatment of radiative transfer than SNEC. For our specific comparison, we choose model 
mlSMdot, which was evolved with enhanced mass loss. The set of parameters provided in Dessart et al. (2013) is sufficient to 
generate a MESA progenitor model with similar but not exactly the same as Dessart et al.’s model mlSMdot. The presupernova 
hydrogen-rich envelope mass in our version of model mlSMdot is Mh ~ 8.13 Mq and the model’s radius is i? ^ 789 Rq. 
Dessart et al. (2013) find Mh ~ 7.72 Mq and R ~ 776 Rq. we attribute the differences to our version to the different release 
versions of MESA used by Dessart et al. (2013) and us. Given the differences in the models (and in the codes), we do not expect 
perfect agreement. 

Our SNEC explosions are triggered by a thermal bomb, while Dessart et al. (2013) use a piston. We tried to match the explosion 
parameters used by Dessart et al. (2013) (given in their Table 2) as closely as possible. In particular, we excised the inner 1.5 Mq 
in agreement with the location of the piston in the models of Dessart et al. (2013). We choose the explosion energy in such a 
way that the final total energy of the model is 1.28 x 10®^ erg. To achieve this, we inject 1.44 x 10®^ erg into the innermost 
0.02 Mq over a time of 0.1 sec. In addition, we include 0.081 Mq of ®®Ni, initially uniformly distributed in the mass coordinate 
range 1.5 — 3.5 Mq. We apply boxcar smoothing, as for the other models from our study. Although we do not know the precise 
distribution of ®®Ni in the model of Dessart et al. (2013), the authors indicate that it was not strongly mixed outwards. We used 
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Figure 17. Comparison of the bolometric light curve of model ml5Mdot of (blue line, Dessart et al. 2013) with the light curve that we obtain with SNEC for a 
similar progenitor model, the same explosion energy, and the same initial ®®Ni mass (black line). The Dessart et al. light curve is generated with their CMFGEN 
non-LTE radiative transfer code that assumes homologous expansion Hillier & Dessart (2012). 


the same prescription for the opacity floor as for all our other models in this study. 

Figure 17 compares our SNEC light curve with the light curve obtained by Dessart et al. (2013) and shown in their Figure 7. 
We find encouraging agreement between the two light curves in the first ~80 days, as well as after day 140 (in the nebular phase). 
The Dessart et al. plateau is 10 — 20 days longer than predicted by SNEC. There are a couple of possible reasons for this difference 
to consider. 

On one hand, the two progenitor models are not identical and the Dessart et al. model has less mass (by ~0.4 Mq) in its 
hydrogen-rich envelope than our model. However, the higher hydrogen-rich mass should result in a somewhat longer plateau 
(e.g. Young 2004; Kasen & Woosley 2009) in our model, which is not what we find in Figure 17. Also, one notes from Figure 17 
that the two light curves would agree much better if the Dessart et al. light curve was shifted back in time by ~5 days. It is not 
clear what would cause such a shift. However, we point out that we find in our SNEC calculations that at day 10 the expansion is 
not yet homologous and more internal energy has yet to be converted into kinetic energy of expansion (this point was previously 
discussed by Dessart & Hillier 2011). Nevertheless, Dessart et al., who assume homologous expansion, map to their radiative 
transfer code already at day ten. The consequences of such an early mapping should be explored. 

On the other hand (Dessart, private communication), the differences may be caused by SNEC’s radiation-hydrodynamics 
treatment, the assumption of LTE, and in the opacities we use (see Section 2). CMFGEN solves the time-dependent radiative 
transfer equations, resolving up to 2000 and more so called super-levels in the frequency domain, which typically represent 
5000 — 10000 atomic levels (see Hillier & Dessart 2012). In optically thin regions, CMFGEN uses line opacities and emissivities 
instead of the Rosseland mean values, used by SNEC. This provides a better description of the effects of iron group line blanketing. 
In order to further isolate potential culprits we have carried out experiments in which we varied the explosion energy, the degree 
of ^®Ni mixing, and the opacity floor. None of these tests produce light curves that are substantially closer to the Dessart et al. 
(2013) light curves than what is shown in Figure 17. 




